function r=earthRadius(latitude)
theta = toRadians(latitude);
a=6378.1370*1000 ;
b=6356.7523*1000;
m = radius(power(a, 2) * cos(theta), power(b, 2) * sin(theta));
n= radius(a * cos(theta), b * sin(theta));
r=m/n;
end

